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Abstract 

There has been a recent appreciation of the ecological impacts of zooplanktonic species 
invasions. The North American brine shrimp Artemia franciscana is one such alien invader 
in hyper-saline water ecosystems at a global scale. It has been shown to outcompete native 
Artemia species, leading to their local extinction. We used partial sequences of the 
mitochondrial Cytochrome c Oxidase Subunit 1 (COI or coxl) gene to investigate the genetic 
diversity and phylogeography of A. salina , an extreme halophilic sexual brine shrimp, over 
its known distribution range (Mediterranean Basin and South Africa) and to assess the 
extent of local endemism, the degree of population structure and the potential impact of 
traditional human saltpan management on this species. We also examined the phylogenetic 
relationships in the genus Artemia using COI sequences. Our results show extensive 
regional endemism and indicate an early Pleistocene expansion of A. salina in the Medi¬ 
terranean Basin. Subsequent population isolation in a mosaic of Pleistocene refugia is 
suggested, with two or three refugia located in the Iberian Peninsula. Two instances of long¬ 
distance colonization were also observed. Surprisingly, given its strong phylogeographical 
structure, A. salina showed a signature of correlation between geographical and genetic 
distance. Owing to strong' priority effects', extensive population differentiation is retained, 
despite dispersal via migrant birds and human management of saltpans. The foreseeable 
expansion of A. franciscana is likely to be followed by substantial loss of genetic diversity 
in Mediterranean A. salina. Large genetic divergences between Mediterranean and South 
African A. salina suggest that the latter deserves species status. 
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Introduction 

Lakes and ponds are isolated ecosystems which, in addition 
to the actively dispersed insects and amphibians, include 
a wide array of passively dispersed invertebrates, both 
planktonic and benthic (e.g. copepods, ostracods, rotifers 
and branchiopods), which form important part of the 
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biomass and food webs of these systems. These aquatic 
invertebrates produce diapausing resistant stages in their 
life cycle which allow them to survive 'in situ' during adverse 
periods or drought. There has been a recent realization 
that accidental (e.g. by transport of ballast water of ships) 
or intentional (e.g. through aquaculture, fisheries or pet 
trade) introductions are leading to nonindigenous species 
invasions in aquatic habitats. Invasive events are occurring 
at an alarming rate and often lead to drastic changes in 
ecosystem functioning and species extinctions (Cristescu 
et al. 2001; Bollens et al. 2002; Frisch et al. 2006; Mergeay et al. 
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2006) . The documentation of species invasions amongst 
passively dispersed invertebrates has been made easier by 
the use of polymerase chain reaction (PCR)-based genetic 
analyses, given the abundance of cryptic species amongst 
these organisms (Gomez et al. 2002; Adamowicz & Purvis 
2005; Mergeay et al. 2005; Adamowicz et al. 2007). However, 
conservation and management strategies have largely 
ignored biodiversity losses in passively dispersed aquatic 
invertebrates, which comprise small 'uncharismatic' species 
(Hamer & Brendonck 1997; Belk 1998; Munoz 2007). 

Population structure and phylogeographical analyses 
in passively dispersed continental aquatic invertebrates 
have indicated that, despite their high colonization ability 
via diapausing stages, these organisms often show an 
unexpectedly high degree of genetic differentiation and 
local endemism which reflects a highly reduced level of gene 
flow among established populations (Hebert 1998; Hebert 
et al. 2003b; Penton et al. 2004; De Gelas & De Meester 2005; 
Paland et al. 2005; Ishida & Taylor 2007; see reviews in 
De Meester 1996, 2002). Given this combination of high 
colonization ability and reduced rates of ongoing gene 
flow, two types of phylogeographical scenario have been 
identified in continental zooplankton. In the first, relatively 
old populations, in or around Pleistocene refugia, show a 
deep genetic structure and strong interpopulation diver¬ 
gence leading to a high degree of local endemism (Gomez 
et al. 2000,2007a; Zierold et al. 2007). In the second contrasting 
scenario, little or no geographical structure and reduced 
genetic diversity is found, indicative of recent and rapid 
range expansion into a newly available area, either after 
natural colonization — often postglacial — or following human 
introduction (Weider et al. 1999; Mergeay et al. 2005; Ishida 
& Taylor 2007). A combination of both scenarios can be 
found in the same species depending of the age of its 
habitat across the range (Ishida & Taylor 2007). 

The brine shrimp genus Artemia (Crustacea, Branchiopoda, 
Anostraca) comprises seven recognized sexual species 
and several parthenogenetic lineages with different ploidy 
levels (Gajardo et al. 2002; Baxevanis et al. 2006). Brine 
shrimps are extreme halophilic organisms and can with¬ 
stand a wide salinity range, from 45 g/L to up to 370 g/L 
(Browne & Hoopes 1990), including different types of ionic 
composition (see Van Stappen 2002 for a review). They 
inhabit patchy hypersaline aquatic habitats such as salt 
lakes and saltworks (salines or saltpans) worldwide in 
both inland and coastal localities, except Antarctica (Trian- 
taphyllidis et al. 1998). These anostracans also constitute 
an important food source for aquatic birds (Sanchez et al. 

2007) . Most Artemia strains can reproduce both by ovovivi- 
parity (producing immediately hatching eggs) and oviparity 
(producing diapausing cysts). The reproductive mode is 
often ovoviviparity when environmental conditions are 
favourable. Under adverse conditions, they produce resistant, 
diapausing cysts (i.e. encysted embryos) (Criel 2002). These 


resistant forms accumulate in saltpan or lake sediments 
forming diapausing egg banks, allowing the persistence of 
populations during unfavourable periods (Caceres 1997); 
they are also the stage responsible for dispersal via waterfowl 
(Green et al. 2005; Sanchez et al. 2007). 

Artemia salina (Linnaeus, 1758) is the only sexual brine 
shrimp species endemic to the Mediterranean Basin, and it 
is also found in a disjunct area in South Africa. This South 
African population was found to be reproductively com¬ 
patible with Mediterranean isolates, but there were some 
significant morphological and life-history differences between 
both groups (Amat et al. 1995a). In addition, the sexual 
North American species A. franciscana, widely used as live 
food in aquaculture, has recently been introduced to the 
Mediterranean Basin and has quickly replaced native 
Artemia species in many localities (Amat et al. 2005b, 2007). 
Phylogeographical surveys in hypersaline anostracans 
are lacking, but analyses on freshwater sexual anostracan 
species indicate a high degree of endemism with a large 
proportion of genetic diversity being held by local popula¬ 
tions (Ketmaier et al. 2003, 2005). Additionally, A. salina is 
associated with current or past commercially exploited 
saltpans (Barigozzi 1974; Abatzopoulos et al. 2002; Amat 
et al. 2005a).The effect of at least several millennia of 
traditional human management of these ecosystems 
(Kurlansky 2003) on the population genetic structure of this 
organism is currently unknown. Therefore, an analysis of 
the genetic diversity of A. salina would be expected to lead 
to two very different scenarios for the loss of genetic diver¬ 
sity arising from the ongoing A. franciscana invasion: (i) if 
human management was responsible for a recent expansion 
of A. salina in the Mediterranean area during historical times, 
or had significantly impacted population structure in this 
species, then the biodiversity loss associated to the expansion 
of A. franciscana would be reduced; and (ii) in contrast, if 
A. salina still holds substantial endemism in its local 
populations despite traditional human management, then 
A. franciscana invasion, together with habitat loss, could 
entail a rapid and significant loss of genetic diversity in this 
organism. For these reasons, an assessment of the genetic 
divergence, population structure and local endemism in 
A. salina is urgently required. 

Here, we describe the geographical distribution of 
mitochondrial genetic diversity of the sexual brine shrimp 
A. salina in order to infer its evolutionary and biogeographical 
history (including the potential effects of human traditional 
management of saltpans). We investigate sequence variation 
using phylogenetic and phylogeographical methods in a 
large set of samples covering the entire native distribution 
of this anostracan. Our results provide one of the most 
complete surveys on the phylogeography of an anostracan, 
showing the high level of endemism and strong genetic 
structure in their isolated populations. This work also 
provides important information for management policies 
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Table 1 Populations of A. salina sampled in this study; n, number of samples; n (standard deviation), nucleotide diversity; H (standard 
deviation), gene diversity; *, haplotypes shared between Donana area and Ebro Delta; +, haplotypes shared between Tunisia and Libya; 
1, haplotypes shared between Portuguesas saltpan and Sardinian saltpans ; 1 , collectively known as the Santa Pola saltpans. 2 , collectively 
known as the Sanlucar saltpans. In the column 'Source' we indicate the source of the samples as CB, samples obtained from the cyst bank 
of the Instituto de Acuicultura de Torre de la Sal; AD, adults sampled in saltpans; BE, cysts collected from bird excreta 


ID {n) 

Site 

Latitude 

Longitude 

71 

H 

No. of haplotypes 
(Hap. codes) 

Source 

EBR (10) 

Spain: Tarragona: Ebro Delta, 

La Trinitat salterns 

40°34'N 

00°41'E 

0.0018 (0.0014) 

0.71 (0.11) 

4 (AS05*,27*,31,32) 

BE 

CAM (8) 

Spain: Mallorca: Salobrar de 

Campos salterns 

39°21'N 

03°00'E 

0.0004 (0.0006) 

0.25 (0.18) 

2 (AS39,40) 

CB 

BRAS (13) 

Spain: Alicante: Bras del 

Port 1 salterns 

38°11'N 

00°36'W 

0.0034 (0.0022) 

0.89 (0.07) 

8 (AS01,04,09,10, 
11,12,13,14) 

AD 

BON (13) 

Spain: Alicante: Bonmati 1 salterns 

38°10'N 

00°37'W 

0.0023 (0.0016) 

0.65 (0.11) 

4 (AS01,04,25,26) 

BE 

MAT (14) 

Spain: Alicante: La Mata Lagoon 

38°02'N 

00°41'W 

0.0011 (0.0010) 

0.39 (0.16) 

4 (AS01,02,03,04) 

BE 

ROS (7) 

Spain: Cordoba: Los Rosales salterns 

37°53'N 

04°46'W 

0.0000 (0.0000) 

0.00 (0.00) 

1 (AS33) 

CB 

SCA (7) 

Spain: Jaen: San Carlos salterns 

37°52'N 

03°40'W 

0.0000 (0.0000) 

0.00 (0.00) 

1 (AS34) 

CB 

PIN (13) 

Spain: Murcia: San Pedro del 

Pinatar salterns 

37°49'N 

00°45'W 

0.0039 (0.0025) 

0.79 (0.11) 

7 (AS01,15,16,17, 

18,19,20) 

AD 

MAL (7) 

Spain: Granada: La Malaha salterns 

37°06'N 

03°43'W 

0.0014 (0.0012) 

0.71 (0.18) 

4 (AS35,36,37,38) 

CB 

DON (13) 

Spain: Huelva: Donana National 

Park salterns 

36°52'N 

06°21'W 

0.0059 (0.0036) 

0.75 (0.07) 

4 (AS05*,06,07,08) 

AD 

POR (12) 

Spain: Huelva: Portuguesas 2 salterns 

36°52'N 

06°20'W 

0.0118 (0.0067) 

0.88 (0.06) 

6 (AS06,07,27*,28,291,30) 

AD 

ROC (11) 

Spain: Huelva: El Rocio 2 salterns 

36°48'N 

06°20'W 

0.0035 (0.0023) 

0.61 (0.10) 

3 (AS06,07,27*) 

BE 

CER (16) 

Spain: Almeria: Cerrillos salterns 

36°41'N 

02°40'W 

0.0010 (0.0009) 

0.35 (0.15) 

4 (AS21,22,23,24) 

AD 

SAH (8) 

Tunisia: Sahline: COTUSAL salterns 

35°45'N 

10°45'E 

0.0029 (0.0021) 

0.86 (0.11) 

5 (AS47,52+,53+,54,55) 

CB 

ADH (8) 

Tunisia: Sebkha El Adhibet: 

Mines Maghreb salterns 

33°06'N 

11°19'E 

0.0039 (0.0027) 

0.96 (0.08) 

7 (AS45,46+,47, 
48,49,50,51) 

CB 

JAD (8) 

Morocco: Oualidia: Souzama salterns 

32°53'N 

08°50'W 

0.0016 (0.0013) 

0.75 (0.14) 

4 (AS41,42,43,44) 

CB 

GAR (8) 

Algeria: Ain Beida: Garaet et 

Tarf salt lake 

35°40'N 

07°09'E 

0.0017 (0.0014) 

0.71 (0.12) 

3 (AS56,57,58) 

CB 

WAD (9) 

Egypt: Wadi el Natrun salt lake 

30°24'N 

30°18'E 

0.0006 (0.0007) 

0.39 (0.16) 

2 (AS59,60) 

CB 

LIB (8) 

Libya: Lake Mandara 

26°41'N 

13°18'E 

0.0023 (0.0018) 

0.75 (0.14) 

4 (AS46+,52+,53+,61) 

CB 

CYP (8) 

Cyprus: Larnaca salt lake 

34°53'N 

33°37'E 

0.0000 (0.0000) 

0.00 (0.00) 

1 (AS62) 

CB 

SGI (7) 

Italy: Sardinia: Sta. Gilla salterns 

39°12'N 

08°59'E 

0.0107 (0.0066) 

0.80 (0.13) 

4 (AS291,63,64,65) 

CB 

MOL (8) 

Italy: Sardinia: Molentargius lagoon 

39°12'N 

09°02'E 

0.0091 (0.0056) 

0.86 (0.11) 

5 (AS291,64,66,67,68) 

CB 

MES (8) 

Italy: Sicily: Maria Estela salterns 

37°59'N 

12°31'E 

0.0004 (0.0006) 

0.25 (0.18) 

2 (AS69,70) 

CB 

VEL (8) 

South Africa: Saldanha Bay: 

Veldrif salterns 

33°22'S 

19°02'E 

0.0091 (0.0056) 

0.64 (0.18) 

4 (AS71,72,73,74) 

CB 


with wider implications for the conservation biology of 
hypersaline habitats in the Mediterranean Basin. 

Materials and methods 

Samples and study area 

We sampled 24 localities where Artemia salina had been 
previously described, covering its whole geographical 
distribution (Table 1). Most samples were obtained as dry 
cysts from the 'cyst-bank' of the Instituto de Acuicultura de 
Torre de la Sal (CSIC, Castellon, Spain) (see Table 1 for 
details). Additional samples were from adults collected in 
the field or cysts extracted from bird excreta, hatched and 
reared to adulthood in the laboratory as part of a morpho¬ 
logical study (Table 1). Cysts and adults were preserved in 


100% ethanol until needed for genetic study. We sequenced 
7-16 individuals from each population. 

DNA extraction, PCR and sequencing 

For DNA isolation of adult specimens we used modified 
digestion temperature (56 °C) and RNAse addition (25 mg/ 
mL) for the Cetyl Trimethyl Ammonium Bromide (CTAB) 
protocol published by Palumbi (1996) and Bossier et al. 
(2004). For DNA isolation of cysts, we used an alkaline lysis 
protocol optimized for zooplanktonic diapausing eggs 
(Montero-Pau et al. in press). 

We used specific Artemia primers designed in the same 
position as primers LCO1490/HCO2198 (Folmer et al. 1994) 
to amplify a fragment of the mitochondrial Cytochrome 
c Oxidase Subunit I (COI) gene. These new primers, 1/ 
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2COI_Fol-F (5-ATT CTA CG A ATC AC A AGG ATA TTG G-3') 
and 1 /2COI_Fol-R (5'-TAC ACT TCA GGA TGG CCA AAA 
AAT CA-3'), provided successful amplification and good 
quality sequences. A fragment of 710 base pairs (bp) was 
amplified by PCR under the following conditions: a cycle 
of 3 min at 94 °C, followed by 35 cycles of 45 s at 94 °C, 60 s 
at 55 °C and 60 s at 72 °C, with a final step of 5 min at 72 °C. 
PCR products were purified and subsequently sequenced 
using the BigDye Terminator Sequencing Ready Reaction 
version 3.1 kit (Applied Biosystems) following manufacturer 
instructions, on an ABI 3130x1 and ABI 3730x1 automated 
sequencer. Nucleotide sequences were edited by hand 
and aligned using sequencher™ version 4.5 (Gene Codes 
Corp., © 1991-2005) and further revised by eye. All sequences 
were deposited in GenBank (Accession Nos DQ426827- 
DQ426858 and EU543444-EU543485). 

Genetic analyses 

The A. salina sequences obtained were collapsed into haplo- 
types using collapse version 1.2 prior to genetic analyses 
(Posada 2005). For interspecific phylogenetic analyses, 
we downloaded all available Artemia COI sequences from 
GenBank to add to our A. salina dataset with three different 
objectives: (i) to assess whether our A. salina samples formed 
a monophyletic unit; (ii) to investigate the phylogenetic 
relationships of the geographically divergent A. salina 
population from South Africa (Veldrif, VET); and (iii) to 
explore the usefulness of the COI gene for DNA barcoding 
purposes in Artemia. We used ModelGenerator version 
0.63 (Keane et al. 2006) to select the best model of sequence 
evolution and used the model and estimated parameters 
in phyml version 2.4.4 (Guindon & Gascuel 2003) for 
phylogenetic reconstructions by Maximum Fikelihood (ME). 
Neighbour-joining (NJ) trees were reconstructed in mega 
version 4 (Tamura et al. 2007), using evolutionary distances 
computed with the Maximum Composite Fikelihood 
method using the Tamura-Nei substitution model (Tamura 
et al. 2007). The robustness of the branches for both ME and 
NJ analyses was assessed with 1000 bootstrap replications. 

Similarly, the phylogenetic relationships between Medi¬ 
terranean A. salina haplotypes were reconstructed using 
the ME approach with phyml incorporating the parameters 
and evolutionary model selected with ModelGenerator. 
In addition, we performed a Nested Clade Phylogeographical 
Analysis (NCPA) of haplotype data (see Templeton 2007). 
We constructed a haplotype network using tcs version 
1.21 (Clement et al. 2000), which follows the statistical 
parsimony algorithm described in Templeton et al. (1992). 
Next, we built a nested design following the guidelines 
provided by Templeton et al. (1995). Finally, we tested for 
the existence of geographical associations of the haplotype 
clades in the nested design using GeoDis version 2.5 
(1 000 000 permutations; Posada et al. 2000). The evolutionary 


and phylogeographical patterns responsible for the geogra¬ 
phical distribution of the clades giving significant results 
were identified following the inference key (updated 11th 
November 2005) provided by these authors with GeoDis. 
The haplotypes from the South African population were 
excluded from the NCPA as they showed more than 40 
substitutions with respect to Mediterranean haplotypes, 
well above the 95% confidence criterion. 

Average genetic distances between Artemia taxa, and 
between the A. salina haplotype higher level clades, obtained 
in the nested design, were computed by mega version 4.0 
(Tamura et al. 2007), using uncorrected p-distances and 
Kimura-2 parameter distances (K-2P) in order to compare 
our results with previous analyses reported in another 
zooplanktonic organisms. Estimates of population genetic 
diversity in A. salina (n, nucleotide diversity; H, gene diversity, 
(]) ST values, genetic variability distribution, and haplotype 
frequencies) were performed using arlequin 2.000 
(Schneider et al. 2000). 

Isolation-by-distance (IBD) analysis between A. salina 
populations was carried out using ibdws version 3.14 
(Jensen et al. 2005), which performs a Mantel test and carries 
out Reduced Major Axis (RMA) regression analysis. To 
visualize the relationship between genetic and geo¬ 
graphical distances, pairwise population c)) ST values were 
plotted against geographical distance (Log geographical 
distance). 

Results 

Phylogenetic relationships in Artemia 

After excluding 15 published Artemia sequences downloaded 
from the GenBank database, which were 561 bp long, and 
trimming longer sequences we obtained a 603 bp long 
alignment for all 98 COI sequences with no indels (alignment 
available from the authors). The evolutionary model selected 
by ModelGenerator was TrN + I + G (I = 0.62; G = 4.20) 
with a transition/transversion ratio = 5.76 and base freq¬ 
uencies A = 20.06%, C = 24.95%, G = 20.96% and T = 34.02%. 
Both methods of phylogenetic reconstruction (ML and NJ) 
resulted in similar topologies and supported the monophyly 
of all species, including all of our A. salina haplotypes, with 
minimum support values of 73% and 91%, depending on 
the method of analysis used (Fig. 2). Support was also 
found for a monophyletic cluster formed by A. tibetiana, 
A. urmiana, an undescribed bisexual population from 
Kazakhstan, and a diploid parthenogenetic lineage, with 
A. urmiana as a sister taxa to the parthenogenetic lineage 
and the bisexual population from Kazakhstan. Strong 
support was found for a sister status of Mediterranean 
A. salina and the group formed by the sequences from the 
South African population (VEL). However, South African 
and Mediterranean lineages were highly divergent, with 
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Table 2 Genetic pairwise distances of the COI mitochondrial gene between: (a) seven Artemia taxa and the Artemia population from Veldrif 
(VEL); (b) high-level clades of A. salina of the NCPA. The lower matrix shows the p-distances values and the upper matrix shows the K-2P 
correction values 

(a) 



A. 'Veldrif' 

A. franciscana 

A. tibetiana 

A. urmiana 

A. sinica 

A. persimilis 

A. salina 

Diploid parth. 

A. 'Veldrif' 

_ 

0.225 

0.255 

0.243 

0.226 

0.265 

0.115 

0.240 

A. franciscana 

0.189 

— 

0.185 

0.192 

0.184 

0.227 

0.241 

0.200 

A. tibetiana 

0.209 

0.160 

— 

0.041 

0.180 

0.238 

0.248 

0.044 

A. urmiana 

0.201 

0.165 

0.040 

— 

0.190 

0.241 

0.237 

0.021 

A. sinica 

0.188 

0.159 

0.155 

0.162 

— 

0.187 

0.247 

0.183 

A. persimilis 

0.218 

0.193 

0.199 

0.201 

0.163 

— 

0.248 

0.228 

A. salina 

0.104 

0.199 

0.203 

0.197 

0.201 

0.206 

— 

0.247 

Diploid parth. 

0.199 

0.171 

0.043 

0.021 

0.156 

0.192 

0.202 

— 

(b) 


3-1 

3-2 

3-3 

3-4 

3-5 

3-6 

3-7 


3-1 

_ 

0.026 

0.025 

0.027 

0.017 

0.022 

0.023 


3-2 

0.025 

— 

0.012 

0.018 

0.018 

0.015 

0.012 


3-3 

0.025 

0.012 

— 

0.016 

0.018 

0.013 

0.011 


3-4 

0.027 

0.018 

0.016 

— 

0.023 

0.018 

0.017 


3-5 

0.017 

0.018 

0.017 

0.022 

— 

0.014 

0.016 


3-6 

0.022 

0.015 

0.013 

0.018 

0.014 

— 

0.013 


3-7 

0.023 

0.012 

0.011 

0.016 

0.016 

0.013 

— 



distances of 10.4% (p-distance; see Table 2a). The haplotypes 
from Mexico and Cuba (which were named A. sp. in their 
GenBank record) form a monophyletic grouping within the 
diversity of A. franciscana. 

Pairwise genetic distances between the seven Artemia 
species and the South African population (VEL) of A. salina 
are shown on Table 2a. The p-distances ranged from 2.1%, 
between the diploid parthenogenetic strain and A. urmiana, to 
a 21.8%, between A. persimilis and A. salina from South Africa. 

Haplotype and nucleotide diversity in A. salina 

The length of the A. salina alignment was 617 bp, with the 
232 individuals sequenced yielding 74 haplotypes. No 
insertions, deletions, stop codon or ambiguities were present. 
There were 115 variable sites, 88 of them parsimony infor¬ 
mative. The alignment had a total of 107 synonymous and 
nine nonsynonymous substitutions. 

Most haplotypes were found within single populations, 
with 12 haplotypes shared among two or more populations. 
The average for each population was 3.87 haplotypes, with 
9.67 individuals sequenced on average for each population 
(see Table 1). The ROS, SC A and CYP populations were 
found to be fixed for a single haplotype. The relative fre¬ 
quencies and accession number of each haplotype are listed 
in Appendix I. 

The values of n and H, as well as details of the haplotypes 
present in each population, are given in Table 1. The 7t-values 


ranged from 0.0000 to 0.0118. The highest value was found 
in POR, but SGI, MOL and VEL also showed high 71-values 
compared with the rest of the populations. In general, 
H-values were high, with a maximum in 0.96 and an average 
of 0.58. 

Most pairwise ()) ST comparisons between A. salina popu¬ 
lations were highly significant (Table 3). Non-significant 
(|) ST values involved two groups of populations within a 
reduced geographical area (Table 3 and Fig. 1), one formed 
by MAT, BRAS and BON in eastern Spain, and another 
formed by DON, ROC and POR in the Donana area of 
southern Spain. A third group involved the two Tunisian 
(SAH and ADH) and the Libyan (LIB) populations. 
Additionally, several pairs of distant populations showed 
nonsignificant (]) ST : two southern Spanish localities (DON 
and POR) were not significantly different from EBR (northern 
Spain), a locality more than 700 km apart where A. salina 
has only been identified recently (see discussion); and one 
of these southern Spanish populations (POR) also showed 
nonsignificant (]) ST values when compared to the two 
Sardinian populations (SGI and MOL). 

Mitochondrial phylogeography of Mediterranean 
A. salina 

We found a total of 70 Mediterranean A. salina haplotypes 
with a maximum divergence between haplotypes of 3.4% 
and 4.5%, p-distance and K-2P, respectively (results not 
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Table 3 Pairwise (J) ST values obtained with a Tamura & Nei evolutionary model (similar values were obtained with K-2P model). Values in bold indicate NO statistical significance (P > 0.01) 



MAT 

DON 

BRAS 

PIN 

CER 

BON 

ROC 

POR 

EBR 

ROS 

SCA 

MAL 

CAM 

JAD 

ADH 

SAH 

GAR 

WAD 

LIB 

CYP 

SGI 

MOL 

MES 

DON 

0.795 























BRAS 

-0.005 

0.741 






















PIN 

0.094 

0.723 

0.077 





















CER 

0.929 

0.875 

0.871 

0.868 




















BON 

0.031 

0.761 

-0.028 

0.112 

0.897 



















ROC 

0.875 

-0.004 

0.816 

0.798 

0.924 

0.840 


















POR 

0.699 

0.049 

0.653 

0.637 

0.798 

0.667 

0.164 

















EBR 

0.899 

0.096 

0.822 

0.803 

0.942 

0.855 

0.324 

0.166 
















ROS 

0.781 

0.773 

0.541 

0.557 

0.953 

0.621 

0.876 

0.644 

0.921 















SCA 

0.593 

0.744 

0.286 

0.379 

0.945 

0.346 

0.860 

0.607 

0.909 

1.000 














MAL 

0.795 

0.637 

0.620 

0.638 

0.916 

0.688 

0.781 

0.494 

0.814 

0.878 

0.828 













CAM 

0.841 

0.804 

0.662 

0.647 

0.952 

0.731 

0.890 

0.691 

0.927 

0.958 

0.938 

0.892 












JAD 

0.940 

0.785 

0.882 

0.868 

0.966 

0.906 

0.871 

0.659 

0.901 

0.961 

0.957 

0.926 

0.960 











ADH 

0.869 

0.526 

0.793 

0.777 

0.925 

0.821 

0.681 

0.363 

0.673 

0.873 

0.856 

0.761 

0.886 

0.821 










SAH 

0.900 

0.612 

0.827 

0.813 

0.942 

0.854 

0.751 

0.450 

0.762 

0.914 

0.902 

0.827 

0.919 

0.847 

0.096 









GAR 

0.934 

0.763 

0.872 

0.856 

0.958 

0.899 

0.858 

0.620 

0.888 

0.955 

0.950 

0.915 

0.954 

0.924 

0.828 

0.871 








WAD 

0.957 

0.734 

0.903 

0.889 

0.972 

0.927 

0.840 

0.589 

0.900 

0.983 

0.981 

0.941 

0.978 

0.949 

0.830 

0.882 

0.949 







LIB 

0.914 

0.646 

0.842 

0.829 

0.950 

0.870 

0.780 

0.487 

0.800 

0.934 

0.925 

0.859 

0.936 

0.877 

0.189 

0.001 

0.891 

0.906 






CYP 

0.948 

0.713 

0.865 

0.842 

0.970 

0.903 

0.850 

0.526 

0.900 

1.000 

1.000 

0.950 

0.988 

0.957 

0.780 

0.851 

0.954 

0.969 

0.894 





SGI 

0.813 

0.473 

0.756 

0.739 

0.868 

0.777 

0.608 

0.158 

0.588 

0.772 

0.749 

0.664 

0.801 

0.751 

0.511 

0.586 

0.715 

0.746 

0.625 

0.662 




MOL 

0.775 

0.179 

0.711 

0.691 

0.853 

0.734 

0.368 

-0.028 

0.317 

0.735 

0.703 

0.584 

0.772 

0.731 

0.416 

0.518 

0.687 

0.693 

0.565 

0.626 

0.088 



MES 

0.963 

0.782 

0.911 

0.898 

0.972 

0.935 

0.887 

0.524 

0.924 

0.991 

0.990 

0.953 

0.984 

0.958 

0.857 

0.899 

0.948 

0.975 

0.922 

0.986 

0.405 

0.582 


VEL 

0.985 

0.972 

0.979 

0.978 

0.986 

0.982 

0.977 

0.958 

0.979 

0.981 

0.981 

0.977 

0.981 

0.976 

0.973 

0.975 

0.978 

0.981 

0.977 

0.981 

0.962 

0.963 

0.980 
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Fig. 1 Map of the Mediterranean Basin showing the location of the 23 Mediterranean Artemia salina populations sampled (see Table 1 for 
population codes and details). Note that POR, DON, ROC and BON, BRAS, MAT are different populations, but due to the scale of the figure 
they appear as single points. The figure also shows the phylogenetic relationships between Mediterranean A. salina COI haplotypes overlaid 
on their geographical distribution and potential Pleistocene refugia (areas joining sampling locations with same colour coding as clades). 
Thin lines indicate long-distance migration events identified in our analysis. The haplotypes marked with asterisks are found in these 
distant areas. The tree is a midpoint rooted maximum likelihood (ML) tree (see text for details) with branch support obtained by 1000 
bootstrap pseudo-replicates (numbers next to branches). Only selected bootstrap values over 70% are shown. 


shown). The haplotype network provided by tcs was used 
to construct a nested design (Fig. 3), which was structured 
into 22 first-level clades, 14 second-level clades and seven 
third-level clades. In addition, three second-level clades 
involved unsampled ancestral haplotypes. Out of the 26 
clades with geographical and genetic information used in 
GeoDis, 13 showed statistically significant geographical - 


genetic associations, including the total cladogram 
(Table 4). The inference for the total cladogram is a global — 
Circummediterranean — range expansion with subsequent 
long-distance colonization and fragmentation for clade 3-2 
(the only one in which the Dn value significantly reversed 
from the Dc value) which involves the two most Eastern 
populations (Egypt and Cyprus). Three clades (1-8,2-8 and 
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A. Satina 


A . sp. Kazah start' 


A. sp. Veldrif 
A. p&rsimilis 


A. paithenogen&tica 
A. tibetia 


A. sp. 'Mexico, Cuba 


A. frandscana 

A sinlCB 


0.05 


Fig. 2 Phylogram of the genus Artemia using Maximum 
Likelihood (ML, below branches) and neighbour-joining (NJ, 
above branches) methods. See text for details. Selected bootstrap 
support values over 70% are shown. Undetermined Artemia taxa 
are indicated as A. sp. followed by their geographical origin as 
reported in the GenBank database. 


3-7) indicate more recent episodes of range expansion. All 
of them involve Spanish populations and four haplotypes 
from the Sardinian populations (AS63, AS64, AS66 and 
AS67). In three other clades (2-5, 3-2 and 3-6), allopatric 
fragmentation was inferred. They are formed by unique 
private haplotypes found in single populations either from 
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islands (Mallorca and Cyprus) and the mainland, or clades 
separated by the Mediterranean Sea (Spain, Morocco and 
Egypt). Restricted gene flow with IBD was inferred in 
two clades (2-6, and 3-5) formed by populations more than 
300 km apart (DON and EBR, or MAT and SC A). Finally, 
clade 2-1 showed evidence of long-distance colonization 
involving populations from southwest Spain (Donana area, 
DON; ROC; POR) and Italian islands (Sardinia, MOL; SGI; 
Sicily, MES). 

The phylogenetic and geographical information is repre¬ 
sented in Fig. 1, identifying the geographical distribution 
of the main clades in the Mediterranean and those geo¬ 
graphical areas with unique haplotypes obtained with 
our data set. Only the most divergent lineages show high 
phylogenetic bootstrap values (3rd level clades 3-1 and 3-2, 
as well as haplotypes from populations JAD, GAR and 
clade 2-1), as is usual in intraspecific phylogenetic trees. 

The average pairwise distances values between haplotypes 
from third level Mediterranean clades ranged between 1.1 
and 2.7% for both p-distance and K-2P (see Table 2b). The 
maximum value was found between clades 3-1 (the CER 
population from SE Spain) and 3-4 [involving the Algerian, 
Sicilian, Sardinian and southern Spain (DON, ROC, POR) 
populations]. 

IBD tested in the Mediterranean haplotype data set showed 
two groups of points, but a significant positive correlation 
of genetic and geographical distances (Z = 560.7279, r = 0.50, 
one-sided P < 0.0001 from 30 000 randomizations) and 
a value of R 2 = 0.246 for RMA regression analysis were 
obtained (Fig. 4). 

Discussion 

Phylogeography and local endemism of Artemia salina in 
the Mediterranean Basin 

Artemia salina populations hold substantial haplotype diver¬ 
sity and are highly structured genetically with an important 
component of regional and local endemism. Our results 
strongly suggest that gene flow is not an important factor 
between A. salina populations. Population structure, instead, 
reflects a phylogeographical pattern of persistent founder 
events ('priority effects'), as found in many other passively 
dispersing continental aquatic invertebrates (see De Meester 
et al. 2002). Overlaying this strong genetic structure, we found 
some evidence of range expansions and a few long-distance 
dispersal events, with possible cases of colonization by immi¬ 
grants derived from populations hundreds of km apart. 

To obtain an approximate time frame for the oldest 
expansion event supporting the inference by the NCPA 
(i.e. the initial range expansion of A. salina in the Mediter¬ 
ranean), we estimated divergence times using the reported 
COI molecular clock for snapping shrimps (1.4% sequence 
divergence per million year; Knowlton & Weigt 1998). 
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Fig. 3 Statistical parsimony network showing the nested design used in NCPA. Only second (continuous thin line), third (continuous thick line) and statistically informative lstlevel clades 
(discontinuous thin line) are shown. Sizes of ellipses are proportional to the number of individuals included (ellipses with numbers bottom left). The populations where each haplotype 
was found are shown next to the haplotype code. Lines linking haplotypes indicate single substitutions. Small black circles indicate missing haplotypes. 
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Table 4 NCPA inferences as described in the GeoDis key of 11 November 2005. Different steps and conclusions are indicated in the chain- 
of-inference column. The alternative steps offered by the GeoDis key are included in parentheses. Locality acronyms are explained in Table 1 


Clade 

Chain of inference 

Populations involved 

Demographic event inferred 

1-8 

1,19,20,2,11 YES (12NO) 

ROS,SCA,MAT 

Range Expansion (Contiguous 

Range Expansion) 

1-12 

1,2,11,17NO 

DON,ROC,POR,EBR 

Inconclusive Outcome 

2-1 

1,2,3,5,6,13YES 

MES,SGI,MOL,DON,POR 

Long Distance Colonization possibly 
coupled with subsequent Fragmentation 

2-2 

l,2,ll,17NO 

SAH,ADH,LIB 

Inconclusive Outcome 

2-5 

l,19NO 

ROS,SCA,MAT,CAM 

Allopatric Fragmentation 

2-6 

l,2,3,4NO 

DON,ROC,POR,EBR 

Restricted Gene Flow with Isolation 
by Distance 

2-8 

1,2,11 YES (12,13,14NO) 

MAT,BRAS,BON,PIN 

Range Expansion (Long Distance 
Colonization and/or Past Fragmentation) 

3-2 

l,19NO 

WAD,CYP 

Allopatric Fragmentation 

3-4 

1,19,20,2 INTERIOR STATUS 
CAN NOT BE DETERMINED 

MES,SGI,MOL,DON,POR,GAR 

Inconclusive Outcome 

3-5 

l,2,3,4NO 

MAT,BRAS,BON,PIN,ROS,SCA,CAM 

Restricted Gene Flow with Isolation 
by Distance 

3-6 

1,19,20,2,3,4,9NO 

ADH,MAL,JAD 

Allopatric Fragmentation 

3-7 

1,2,11 YES (12NO) 

DON,ROC,POR,EBR,MOL,SGI 

Range Expansion (Contiguous 

Range Expansion) 

Total Cladogram 

1,2,11YES (12,13YES) 

ALL LOCALITIES 

Range Expansion (Long Distance 
Colonization possibly coupled with 
subsequent Fragmentation) 


Although this calibration assumes a fixed rate of variation 
in shrimps, it is useful enough for our purposes. According 
to this, the split between the seven Mediterranean A. salina 
main clades (see Table 2b and Fig. 3) would have occurred 
during the Early Pleistocene, well before the last glacial 
maximum, between 0.78 and 1.93 million years ago. Gomez 
et al. (2000) found a pattern of regional endemism for the 
salt lake rotifer Brachionus plicatilis on the Iberian Peninsula, 
suggesting two glacial refugia (with an average pairwise 
sequence divergence of 2.8%). Divergences found by De 
Gelas & De Meester (2005) among groups in Europe for the 
freshwater cladoceran Daphnia magna were lower, ranging 
from 0.73% to 1.74%. Similar divergences were found 
by Penton et al. (2004) in North American populations of 
D. obtusa. Divergences in the Mediterranean 3rd-level clades 
of A. salina range from 1.1% to 2.7%, which could be tenta¬ 
tively dated as Pleistocene at both local (Iberian Peninsula) 
and regional (Mediterranean Basin) scales. 

According to (i) the level of genetic differentiation between 
clades; (ii) supported monophyly of the clades; (iii) geogra¬ 
phical isolation; (iv) presence of unique haplotypes; and 
(v) inference of allopatric fragmentation in the NCPA, 
our results suggest that A. salina populations survived 
during several glaciations in at least six suitable areas 
around the Mediterranean Sea (CER, JAD, GAR, Italian 
islands, CYP-WAD and Donana area). The clades found in 
four of these areas show high support in the phylogenetic 



Fig. 4 Data plot of the genetic distance ((J) ST ) vs. geographical 
distance (Log geographical distance in km) of the Mediterranean 
Artemia salina populations, showing the RMA regression line 
(R2 = 0.246). 
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analyses. CYP-WAD is geographically isolated, has high 
support for monophyly respect to the rest of lineages, and 
allopatric fragmentation is inferred by NCPA; and the 
Donana area shows unique haplotypes very differentiated 
from the Italian lineage after inferring long-distance 
colonization coupled with fragmentation by NCPA (see 
Table 4, Figs 1 and 3). The precise location of such long-term 
refugia is difficult to ascertain. Most likely an area, instead 
of a particular saltpan or lake, can be considered as the 
'refugia', given that particular lakes may have a short life 
span. Given this caveat, we propose that at least two refugia 
were located in the Iberian Peninsula (see Figs 1 and 3). 
These areas contain geographically differentiated clades, 
indicating long-term survival within the Iberian Peninsula 
which we propose are located in or around the following 
areas: (i) clade 3-1, southeastern Area (Almeria), which is very 
isolated genetically from the other Iberian populations; 
and (ii) clade 3-7, southwestern region (Donana area). This 
clade displays two significant historical signatures. One 
of them, a range expansion from the Iberian Peninsula 
towards Sardinia, is considered below in detail. The other 
is a signature of restricted gene flow with IBD between the 
Ebro (northeastern) and Donana areas, which share two 
haplotypes. After Donana, the Ebro Delta is the most 
important locality for migratory waterfowl in the Iberian 
Peninsula (Marti & Del Moral 2002). Several direct move¬ 
ments of ringed birds between both areas have been 
documented (Spanish Ringing Office, unpublished data), 
making dispersal of A. salina cysts between them particularly 
likely. Since A. salina has not previously been detected in 
the Ebro Delta population (Amat et al. 1995b), this could 
reflect the establishment of a new population after long¬ 
distance transport from the Donana area. However, the high 
haplotype richness and the presence of unique haplotypes 
in the Ebro Delta population indicate that this population 
could be older, although a larger sample size might reveal 
these haplotypes to be present in the Donana or surrounding 
areas. An additional refugium may possibly have been 
located in the Iberian Peninsula, harbouring the ancestors 
of clade 3-5 — interior and eastern coastal Iberian region 
(Levante Spanish coast and Balearic Islands) — which is 
supported by the high genetic differentiation of clade 3-5 
from the other 3rd level clades. Sub-clades nested in clade 
3-5 (clades 1-8,2-5 and 2-8) have been inferred by NCPA to 
have undergone (i) a 'contiguous range expansion' in clade 
1-8; (ii) 'allopatric fragmentation' between internal popu¬ 
lations (SCA, ROS) and the geographically isolated 
populations MAT and CAM; and (iii) 'range expansion 
with long-distance colonization' for the Eastern Spanish 
coast populations, MAT, BRAS, BON and PIN. These 
inferences are concordant with a colonization route from 
interior saltpans towards the eastern coast of the Iberian 
Peninsula, arriving at the Balearic Islands (see Table 4 and 
Fig. 3 for details). The PIN population is the only Spanish 


eastern coastal population with significant ([) ST values 
compared to the others. The same occurs with the CAM 
population in the Balearics. Hence, the hypothetic colo¬ 
nization route is supported by 'allopatric fragmentation' 
and 'long-distance colonization' when the populations 
CAM and PIN are involved. A pattern of refugia within 
refugia has also been reported in Iberia for other organisms 
(Gomez & Lunt 2007), including salt lake invertebrates 
(Gomez et ah 2000, 2007a), and appears to be a common 
pattern in other Mediterranean peninsulas (Schmitt et ah 
2006; Canestrelli et ah 2007). 

Migratory birds are the main natural dispersal vectors 
of A. salina between patchy habitats (Figuerola et ah 2005; 
Sanchez et ah 2007). When suitable habitats are in expansion, 
therefore, colonization will not necessarily follow a pattern 
involving occupation of adjacent populations (McMaster 
et ah 2007). Previously colonized populations should be very 
resilient to the effects of new migrants due to monopolization 
effects (De Meester et ah 2002), but 'empty habitats' could 
be colonized by distant populations depending on the nature 
of bird movements. We found at least two examples of 
long-distance colonization, and one case of bidirectional 
exchange of migrants between areas. One of these bidi¬ 
rectional events was between the two clades involving 
Sardinia and southwestern Spain (Donana area) populations. 
A relatively old episode of bidirectional migration between 
these areas followed by isolation was inferred from NCPA. 
This genetic structure coincided with a structured geo¬ 
graphical distribution, but 'long-distance colonization' 
and 'range expansion' events were found between Donana 
and Sardinian populations. Surprisingly, two 3 rd -level 
clades were found in both geographical areas (clades 2-1, 
and clade 3-7 shown in Fig. 3). As the network shows, the 
mitochondrial lineage from clade 2-1 (with basal haplotypes 
in Sardinia and related haplotypes in Sicily) may have 
undergone long-distance dispersal and colonization in the 
past from Sardinia towards the Donana area. Similarly, 
the lineage from clade 3-7 (with basal haplotypes in Donana 
and Ebro Delta populations) could have suffered a range 
expansion from the Iberian Peninsula towards Sardinia. 
The presence of new tip haplotypes in each putatively 
colonized area indicates that the colonization event is 
unlikely to be recent (i.e. due to human management in 
the last few millennia). In any case, the colonization from 
Donana to Ebro Delta appears to be more recent than that 
involving Donana and Sardinia. 

Our results rule out the existence of two geographically 
differentiated (eastern and western) A. salina clades in the 
Mediterranean Basin as recently proposed by Baxevanis 
et ah (2006), based on data from a small number of popula¬ 
tions. Our work demonstrates that Mediterranean A. salina 
is subdivided into seven or eight main clades with mostly 
allopatric distributions, reduced gene flow and long-distance 
dispersal with fragmentation. 
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Given the strong genetic structure, the strong priority 
effects and subsequent monopolization likely to influence 
A. salina populations, migration-drift equilibrium as a result 
of restricted gene flow appears unlikely in the Mediterranean 
populations of A. salina. Therefore, the significant correlation 
between genetic and geographical distance (i.e. IBD pattern) 
could instead reflect a pattern of past sequential colonization 
as found in humans (Ramachandran et al. 2005), and also 
suggested to be responsible for the patterns found in several 
marine and freshwater aquatic invertebrates (Gomez et al. 
2007a, b; Mills et al 2007). 

There is an ongoing debate on the usefulness and short¬ 
comings of NCPA (Garrick et al. 2008; Petit 2008; Templeton 
2008). The method has recently attracted criticism (Knowles 
& Maddison 2002; Panchal & Beaumont 2007; Petit 2008). 
In response to these criticisms the robustness of inferences 
using this approach has been further tested (Templeton 
2004), and its use has been defended given its usefulness by 
generating specific testable hypothesis (Garrick et al. 2008). 
NCPA, indeed, offers a unique ability to explore patterns 
relating to complex historical scenarios. The high rate of 
false positives found by Panchal & Beaumont (2007) for 
some of the inferences could have been generated by the 
specific parameters of the simulations tested which are far 
from those found in natural scenarios. We therefore agree 
with Garrick et al. (2008) that NCPA is a good starting 
point To attempt distinguishing processes and events that 
shaped spatial-genetic structuring throughout complex 
evolutionary histories of natural populations', especially 
when cross-validated 'via assessment of multiple independent 
loci, complementary analyses, and/or prior expectations 
in conjunction with complementary analyses and strong 
prior expectations'. One consequence of the application of 
the interpretive key of NCPA to our dataset is noteworthy. 
In Artemia, as in other passively dispersed organisms in 
which migratory birds are the vectors of dispersal, a clear- 
cut distinction between 'contiguous range expansion' and 
'long-distance colonization' is difficult to make, as their 
habitats are naturally isolated and often far from each other. 
So, the inference key provided by Posada and Templeton 
probably does not consider the high potential of dispersal 
and patchy nature of the habitats of these organisms 
(Posada et al. 2002). We suggest that in organisms that inhabit 
patchy and isolated habitats and are passively dispersed, 
the NCPA resulting inferences 'range expansion' and 
'contiguous range expansion' could be considered to 
represent the same phylogeographical and evolutionary 
event at different but continuous scales. For instance, our 
results indicate 'range expansion' at first, second and third 
clade levels, and at the total cladogram level (see results 
section for details). The differences in geographical dis¬ 
tances between the haplotypes involved in these clades are 
large, but the inference is the same. Another case is repre¬ 
sented by an inference of 'contiguous range expansion' in 


clades whose haplotypes are separated by the sea (clade 3-7 
found in southwestern Spain and Italian islands), where 
obviously no populations exist in between. 

Phylogenetic relationships of Artemia 

Inter-specific phylogenetic relationships in Artemia show 
that all the taxa included are monophyletic. In addition, 
a well-supported grouping of A. tibetiana, A. urmiana and 
diploid parthenogens was found (see Fig. 2), in agreement 
with previous work indicating that diploid parthenogenetic 
lineages are closely related to A. urmiana (Baxevanis et al. 
2006). However, a surprising result was the near identical 
sequences of the A. sp. 'Kazakhstan' bisexual isolate (not 
analyzed by Baxevanis et al. 2006) and the diploid par¬ 
thenogens. We suggest that in view of these results, the 
taxonomic status of the central Asian species should be 
re-evaluated. The Mexican and Cuban sequences can be 
considered conspecific with A. franciscana (Tizol-Correa 
2006), suggesting that A. franciscana has also a strong phylo¬ 
geographical structure with a similar depth to that 
found in A. salina , but more exhaustive analyses are 
needed in future studies. Several approaches have been 
proposed to identify different Artemia species, such as 
morphological parameters (Hontoria & Amat 1992; Zhou 
et al. 2003), or mitochondrial DNA restriction fragment 
length polymorphisms (Bossier et al. 2004). But, as has 
been suggested for most animal groups (Hebert et al. 2003a), 
our results support the use of COI as a valuable 'DNA 
bar-coding' tool in Artemia. Given the high interspecific 
diversity, it could also be used to assess contamination of 
Artemia cultures (see Campos-Ramos et al. 2003) using 
high throughput techniques such as TaqMan© Real Time 
PCR already used in other organisms (Fox et al. 2005; 
Walsh et al. 2005). 

A remarkable finding of our work is the large divergence 
between the South African and the Mediterranean A. salina 
sequences. The phylogram shows that the South African 
haplotypes form a well-supported sister lineage to the 
Mediterranean ones, indicating that this population has 
diverged strongly from Mediterranean A. salina populations 
(see Fig. 2). The average sequence divergence between 
populations of both areas (see Table 2a) is within the range 
of the interspecific divergences found in COI for crustaceans 
(Costa et al. 2007), and it is much larger than the divergence 
between other sexual established Artemia species such as 
A. tibetiana, A. sinica, A. urmiana and the diploid parthenogens 
(Baxevanis et al. 2006). Although the VET population showed 
reproductive compatibility with a Mediterranean popula¬ 
tion, clear morphological and life-history differences could 
be identified between South African and Mediterranean 
isolates (Amat et al. 1995a). Reproductive compatibilities 
have also been reported between other sexual Artemia 
species such as A. sinica, A. urmiana and A. sp. 'Kazakhstan' 
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(Triantaphyllidis et al. 1998). Therefore, the taxonomic status 
of the South African A. salina should be re-evaluated. 

Consequences for conservation and population 
management 

Saltpans are known to have existed and been managed by 
humans for at least 2000 years in the Mediterranean for 
evaporative salt production (Consejeria de Medio Ambiente 
2004), although others are far more recent. Our results on 
the structure and the age of the clades in A. salina suggest 
that traditional human management of saltpans has not 
played a significant role in the genetic structure of its 
populations, which we would have expected to result in 
strong signals of recent expansions or gene flow. Even in 
those cases of long-distance dispersal and colonization, the 
high haplotype diversity and number of unique haplotypes 
appears to rule out recent episodes of anthropogenic-derived 
dispersal. In addition, and given the level of differentiation 
and genetic divergence between geographical isolates, our 
results indicate that A. salina was either pre-existent in the 
natural saltpans or lagoons that later became saltworks, or 
colonized these from nearby natural populations. Therefore, 
traditional forms of saltpan management seem to have 
aided in the maintenance of the genetic diversity of A. salina 
populations; or at least, not to have resulted in severe losses 
of genetic diversity. In contrast to this low impact of 
traditional salt management, many saltpans are currently 
being abandoned or transformed into intensive aquaculture 
projects due to decreased economic returns of salt production 
(Kurlansky 2003; Amat et al. 2005a). After salt production 
is discontinued, saltpans usually become unsuitable for 
Artemia eventually, owing to conversion for aquaculture 
facilities or tourist development (Alonso-Villalobos et al. 
2001), siltation (Aguilera & Gracia 2004) or due to a 
reduction in salinity due to lack of marine water influx 
or the impact of irrigation in surrounding areas (e.g. at 
CER; M. Paracuellos, personal communication). In addition, 
due to its value as fish food, the North American species 
A. franciscana has been introduced to different areas of 
the Mediterranean and has already colonized saltpans 
from Portugal, Spain, France, Italy and Morocco, causing 
the local extinction of native Artemia species (Amat et al. 
2005a, 2007). 

Our results provide evidence that native Artemia can be 
dispersed over long distances by birds. The invasive A. 
franciscana uses the same dispersal vector (Green et al. 2005). 
Both A. salina and A. franciscana were detected in the Ebro 
Delta population for the first time in 2004 (Amat et al. 2005a), 
whereas A. franciscana appears to have replaced A. salina in 
the main flamingo breeding site in Spain, Fuente de Piedra 
lagoon (Malaga), around the same time (Amat et al. 2007). 
Thus, further studies regarding its dispersal and establish¬ 
ment capacity must be carried out. 


The high diversity of saline environments in the Iberian 
Peninsula is reflected in an exceptional endemic genetic 
diversity as has been documented in studies of salt lake 
fauna, such as rotifers (Gomez et al. 2000,2007a), and aquatic 
coleopterans living in saline streams (Abelian et al. 2007). 
The evidence we provide on the high genetic diversity for 
A. salina in the Iberian Peninsula indicates that there is a 
risk of irreversible loss of genetic diversity caused by the 
extinction of local populations. The presence of the dia- 
pausing eggs or cysts in the Artemia life cycle could be 
exploited to enable the preservation of genetic diversity. 
The establishment of long-term cyst banks for all known 
populations could be an initial step for the conservation of 
A. salina diversity. These cysts would allow the possibility 
to restore lost populations at some time in the future if 
propagule banks are sampled timely. 

Conclusions 

Artemia salina populations hold substantial genetic diversity, 
with a strong phylogeographical structure and high regional 
endemism in the Mediterranean Basin and Iberian Peninsula. 
Our analyses point to an episode of population expansion 
in this area, possibly during the early Pleistocene, followed 
by more localized expansions through recolonization routes. 
These expansions sometimes involved long-distance colo¬ 
nization or isolation and genetic fragmentation. 

Given the strong genetic structure, migration-drift equ¬ 
ilibrium as a result of restricted gene flow in equilibrium 
conditions appears to be unlikely in the Mediterranean 
populations of A. salina. The signature of correlation between 
genetic and geographical distance could instead reflect a 
pattern of sequential colonization due to an ancient 
founder effect. 

Our results indicate that the A. salina population from 
South Africa is highly differentiated genetically from the 
Mediterranean populations and might deserve specific status. 

Finally, our results indicate that traditional human saltpan 
management used to be compatible with the maintenance 
of population genetic diversity in this organism. However, 
recent changes of saltpan use, such as aquaculture, which 
promotes the introduction of the alien species A. franciscana, 
as well as habitat loss, pose serious threats for the genetic 
diversity of A. salina populations in the Mediterranean area. 
Due to these threats, A. salina requires urgent measures to 
conserve its remaining genetic diversity throughout its range. 
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Appendix I Summary of the relative frequencies of haplotypes, number of individuals per population and GenBank Accesion no 


A. salina populations 


MAT DON BRAS PIN CER BON ROC POR EBR ROS SCA MAL CAM JAD ADH SAH GAR WAD LIB CYP SGI MOL MES VEL GenBank 

Haplotype (n = 14) (n * 13) (n = 13) (n = 13) (n = 16) (n = 13) (n = 11) (n = 12) (n = 10) (n = 7) (n = 7) (n = 7) (n = 8) (n = 8) (n = 8) (n = 8) (n = 8) (n = 9) (n = 8) (n = 8) (n = 7) (n = 8) (n = 8) (n = 8) Acc.Num. 


AS01 0.7860 - 0.3080 0.4620 - 0.5380 

AS02 0.0714 ----- 

AS03 0.0714 ----- 

AS04 0.0714 - 0.1540 - - 0.3080 

AS05 - 0.2310 - 

AS06 - 0.3080 

A SO7 - 0.3850 - 

AS08 - 0.0769 - 

AS09 - - 0.0760 - 

AS10 - - 0.1540 - 

ASH - - 0.0760 - 

AS12 - - 0.0760 - - - 

AS13 - - 0.0760 - - - 

AS14 - - 0.0760 - - - 

AS15 - - - 0.1540 - - 

AS16 - - - 0.0760 - - 

AS17 - - - 0.0760 - - 

AS10 - - - 0.0760 

ASK) - - - 0.0760 - 

AS20 - - - 0.0760 - - 

AS20 - - - - 0.0620 - 

AS20 - - - - 0.0620 - 

AS20 - - - - 0.8120 - 

AS20 - - - - 0.0620 - 

AS20 _____ 0.0760 

AS20 _____ 0.0760 

AS20 

AS20 - - - - - - 

AS20 - - - 

AS30 

AS30 - - - - 

AS30 - - - - 

AS30 ______ 

AS30 

AS30 

AS30 - - - - - - 

AS30 - - - 

AS30 - - - - - 

AS30 - - 

AS40 - - - - - - 

AS40 

AS40 

AS40 - - - 

AS40 ______ 


0.5000 - 

0.3640 0.2500 - 
0.5450 0.2500 - 


0.0900 0.1670 0.1000 - 

0.0830 - _____ 

0.0830 - 

0.1670 - 

0.3000 - 

0.1000 _____ 

- 1.0000 - 

1.0000 - 

0.5710 - 

_ 0.1430 

0.1 ISO - 

- - - - - 0.1430 - - 

0.8750 
0.1270 - 

- 0.1250 

0.2500 

0.5000 

- - - - 0.1250 


0.4290 0.2500 - 


DQ426827 

DQ426828 

DQ426829 

DQ426830 

DQ426831 

DQ426832 

DQ426833 

DQ426834 

DQ426835 

DQ426836 

DQ426837 

DQ426838 

DQ426839 

DQ426840 

DQ426841 

DQ426842 

DQ426843 

DQ426844 

DQ426845 

DQ426846 

DQ426847 

DQ426848 

DQ426849 

DQ426850 

DQ426851 

DQ426852 

DQ426853 

DQ426854 

DQ426855 

DQ426856 

DQ426857 

DQ426858 

EU543444 

EU543445 

EU543446 

EU543447 

EU543448 

EU543449 

EU543450 

EU543451 

EU543452 

EU543453 

EU543454 

EU543455 
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Appendix I Continued 


A. salina populations 

MAT DON BRAS PIN CER BON ROC POR EBR ROS SCA MAL CAM JAD ADH SAH GAR WAD LIB CYP SGI MOL MES VEL GenBank 

Haplotype (n = 14) (n = 13) (n = 13) (n = 13) (n = 16) (n = 13) (n = 11) (n = 12) (n = 10) (n = 7) (n = 7) (n = 7) (n = 8) (n = 8) (n = 8) (n = 8) (n = 8) (n = 9) (n = 8) (n = 8) (n = 7) (n = 8) (n = 8) (n = 8) Acc.Num. 


AS40 ______________ 0.1250 -------- - EU543456 

AS40 ______________ 0.2500 - - - 0.1250 ----- EU543457 

AS40 ______________ 0.1250 0.2500 ------- - EU543458 

AS40 ______________ 0.1250 -------- - EU543459 

AS40 ______________ 0.1250 -------- - EU543460 

AS50 ______________ 0.1250 -------- - EU543461 

AS50 ______________ 0.1250 -------- - EU543462 

AS50 _______________ 0.3750 - - 0.5000 ----- EU543463 

AS50 _______________ 0.1250 - - 0.2500 ----- EU543464 

AS50 _______________ 0.1250 ------- - EU543465 

AS50 _______________ 0.1250 ------- - EU543466 

AS50 ________________ 0.2500 ------- EU543467 

AS50 ________________ 0.5000 ------- EU543468 

AS50 ________________ 0.2500 ------- EU543469 

AS50 - - - - - - - - - ________ 0.7780 ------ EU543470 

AS60 - - - - - - - - - ________ 0.2220 ------ EU543471 

AS60 - - - - - - - - - _________ 0.1250 ----- EU543472 

AS60 - - - - - - - - - __________ 1.0000 - - - - EU543473 

AS63 - - - - - - - - - ___________ 0.1430 - - - EU543474 

AS64 - - - - - - - - - ___________ 0.1430 0.3750 - - EU543475 

AS65 - - - - - - - - - ___________ 0.2860 - - - EU543476 

AS66 - - - - - - - - - ____________ 0.1250 - - EU543477 

AS67 - - - - - - - - - ____________ 0.1250 - - EU543478 

AS68 - - - - - - - - - ____________ 0.1250 - - EU543479 

AS69 - - - - - - - - - _____________ 0.8750 - EU543480 

AS70 - - - - - - - - - _____________ 0.1250 - EU54348 

AS71 - - - - - - - - - ______________ 0.6250 EU54348 

AS72 - - - - - - - - - ______________ 0.1250 EU54348 

AS73 - - - - - - - - - ______________ 0.1250 EU54348 

AS74 - - - - - - - - - ______________ 0.1250 EU54348 
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